sampleSize=1700;
outputSize=1000000;
lagMax=3; #Max lag for autocor
ro<-seq(length=lagMax, 0, 0); #autocor
Y<-seq(length=outputSize, 0, 0); #output

sampleMean=80.684560000;
sampleSd=0.212428749;

#generate observed samples
#X=rnorm(sampleSize,mean=0, sd=1);

X<-seq(length=sampleSize, 0, 0)
for(i in seq(1, sampleSize, 1))
{
	X[i]=(a[i]-sampleMean)/sampleSd;
}

Ro=acf(X, lag.max=lagMax, type="correlation", plot=FALSE);

#build covariance matrix (inversely), M means N-1
for(i in seq(1, lagMax, 1))
{
	ro[i]=Ro$acf[lagMax+2-i];
}
#print ro
ro;

Ro_N_N=matrix(c(1, ro[1], ro[2], ro[3],  ro[1], 1, ro[1], ro[2], ro[2], ro[1], 1, ro[1], ro[3], ro[2], ro[1], 1), nrow=4, byrow=TRUE);
Ro_1_M=matrix(c(ro[1], ro[2], ro[3]), nrow=1, byrow=TRUE);
Ro_M_1=matrix(c(ro[1], ro[2], ro[3]), nrow=3, byrow=TRUE);
Ro_1_1=matrix(c(1), nrow=1, byrow=TRUE);
Ro_M_M=matrix(c(1, ro[1], ro[2], ro[1], 1, ro[1], ro[2], ro[1], 1), nrow=3, byrow=TRUE);

#initialize output
for(i in seq(1, lagMax, 1))
{
	Y[i]=X[i];
}
mu_matrix=matrix(c(Y[1], Y[2], Y[3]), nrow=3, byrow=TRUE);

temp1=Ro_1_M %*% (solve(Ro_M_M));
temp2=Ro_1_M %*% (solve(Ro_M_M)) %*% Ro_M_1;

#calculate new mu and rio
for(i in seq(4, outputSize, 1))
{
	for(j in seq(1, lagMax, 1))
	{
		mu_matrix[j][1]=Y[i-lagMax-1+j];
	}
	new_mu=temp1 %*% mu_matrix;
	new_ro=1-temp2;
	Y[i]=rnorm(1, mean=new_mu, sd=new_ro);
}